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We study the stability and chaos of three compact objects using post-Newtonian (PN) equations 
of motion derived from the Arnowitt-Deser-Misner-Hamiltonian formulation. We include terms 
up to 2.5 PN order in the orbital part and the leading order in spin corrections. We performed 
numerical simulations of a hierarchical configuration of three compact bodies in which a binary 
system is perturbed by a third, lighter body initially positioned far away from the binary. The 
relative importance of the different PN orders is examined. The basin boundary method and the 
computation of Lyapunov exponent were employed to analyze the stability and chaotic properties 
of the system. The 1 PN terms produced a small but noticeable change in the stability regions of 
the parameters considered. The inclusion of spin or gravitational radiation does not produced a 
significant change with respect to the inclusion of the 1 PN terms. 
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I. INTRODUCTION 

The dynamics of compact objects play an important 
role in the evolution of galaxies and other stellar sys- 
tems. Post-Newtonian (PN) techniques are a useful tool 
for modeling the dynamics of multiple compact objects. 
In astrophysical models, the gravitational radiation is in- 
cluded via an effective force which considers 1 PN and 
2.5 PN corrections, and in some cases stellar dynamical 
friction. Hierarchical three black hole configurations in- 
teracting in a galactic core have been the focus of recent 
studies by several authors. For example, intermediate- 
mass black holes with different mass ratios are consid- 
ered in [IH3] ; simulations of dynamical evolution of triple 
equal-mass supermassive black holes in a galactic nuclei 
were performed in [3|[5]. Other astrophysical applications 
of multiple black hole simulations include three-body 
kicks [6l[7] and binary-binary encounters (see e.g. [8HT2]). 
In the context of the final parsec problem [T31 [2] , three- 
body interactions are considered as a mechanism that can 
drive a binary black holes system to a separation below 
one parsec. 

Since the 1990s there has been an increasing effort to 
detect and study compact objects. The most likely source 
of gravitational waves are binary compact objects. Re- 
cently, it was shown that the probability of more than two 
black holes to interact in the strongly relativistic regime 
is, not surprisingly, very small |15) . For practical pur- 
poses, the creation of gravitational waveform templates 
for gravitational wave detectors is naturally focused on 
binary systems. Binary systems can produced complicate 
waveforms when taking into account spinning black holes 
and eccentric orbits, e.g. [16]. On the other hand, triple 
systems consisting of a binary black hole and a star (e.g., 
a white dwarf) are potential sources of electromagnetic 
counterparts associated with the mergers (T7H2D]. 
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The chaotic behavior of triple systems is well known in 
the Newtonian case (see [UJ and references therein) . For 
binaries, it is known that chaos appears when using cer- 
tain post-Newtonian approximations for spinning binary 
systems (see [22H28] ). In the present work, we studied 
three-body systems with PN methods, where the main 
technical novelty is the inclusion of the 2.5 PN terms in 
the orbital dynamics and leading order of spin-orbit and 
spin-spin terms [29lf3"T] . 

Recently, similar PN techniques were applied to the 
general relativistic three-body problem. Periodic solu- 
tions were studied using the 1 PN and 2 PN approxima- 
tions in [32H33]- Examples of three compact bodies in 
a collinear configuration were considered in [351 , and 
Lagrange's equilateral triangular solution was studied in- 
cluding 1 PN effects in [37 . In [38] , the stability of the 
Lagrangian points in a black hole binary system was stud- 
ied in the test particle limit, where the gravitational radi- 
ation effects were modeled by a drag force. The waveform 
characterization of hierarchical non-spinning three-body 
configurations using up to 2.5 PN terms was presented 

in [32]. 

Close interaction and merger of black holes require nu- 
merical relativistic simulations. The first complete simu- 
lations, using general-relativistic numerical evolutions of 
three black holes, were presented in j4DJ[JT]. These recent 
simulations showed that three compact object dynamics 
display a qualitatively different behavior than Newtonian 
dynamics. In [32], the sensitivity of fully relativistic evo- 
lutions of three and four black holes to changes in the 
initial data was examined, where the examples for three 
black holes are some of the simpler cases already dis- 
cussed in rjOl [IT] . The apparent horizon and the event 
horizon of multiple black holes have been studied in [431 - 
05]. Although fully general-relativistic simulations are 
available, they are constrained by the number of orbits 
and separations between black holes. 

The paper is organized as follows: in Sec. [TTJ we sum- 
marize the equation of motion up to 2.5 post-Newtonian 
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approximation for three spinning bodies. This is followed 
by a discussion of the chaos indicators that we used to 
characterize the triple system. In Sec. |III A[ we describe 
the numerical techniques used to solve the equation of 
motion, and we present some results for test cases. The 
perturbation of a binary system by a third object is pre- 
sented in Sec. |IIIB| where we performed numerical ex- 
periments in order to study the stability and chaos of the 
triple system. The conclusions are presented in Sec. IV 



A. Notation and units 

We employed the following notation: x — (xi) denotes 
a point in the three-dimensional Euclidean space K 3 and 
letters a,b, . . . are the particle labels. We defined f a := 
x- x a , r a := \r a \, h a := f a /r a ; for a ^ b, f ab := x a - 
Xb, r ab := \f ab \ and h ab := f ab /r ab ; here, | • | denotes 
the length of a vector. The mass parameter of the a-th 
particle is denoted by m a with M = J2 a m a- Summation 
runs from 1 to 3. The linear momentum vector is denoted 
by p a . A dot over a symbol, x, means the total time 
derivative, and the partial differentiation with respect to 
x 1 is denoted by <9j. 

In order to simplify the calculations, it is useful to de- 
fine dimensionless variables (see e.g. [46]). We used as ba- 
sis quantities for the Newtonian and post-Newtonian cal- 
culation the gravitational constant G, the speed of light 
c and the total mass of the system M. Using derived 
constants for time r = MG/c 3 , length I — MG/c 2 , lin- 
ear momentum V = Mc, spin S = M 2 G/c and energy 
£ = Mc 2 , we construct dimensionless variables. The 
physical variables are related to the dimensionless vari- 
ables by means of scaling. Denoting with capital letters 
the physical variables with the standard dimensions and 
with lowercase the dimensionless variables. We defined 
for a particle a its position x a :— X a /l, linear momentum 
p a := P a jV and mass m a = M a /M (notice that m a < 1, 
Va). 



II. EVOLUTION METHOD 

A. Equations of motion 

In the ADM post-Newtonian approach, it is possible to 
split the orbital and spin contribution to the Hamiltonian 
(see e.g. [13 US]) 

H(x a ,p a ,S a ) = H(x a ,p a )orb + H(x a ,p a ,S a ) Spin , (1) 

where each component forms a series with coefficients 
which are inverse powers of the speed of light. We in- 
cluded terms up to 2.5 PN contributions where the or- 
bital part for three bodies is given by [3H [3JB SS] 



Here each term of the Hamiltonian is labeled by a super- 
script n that denotes the PN order (powers of c~ 2n ) and 
a subscript which distinguished between the Newtonian 
terms and the PN components. The spin contribution 
forms a power series where the leading order is given by 
[50H53] (see [31] for the next-to-leading order terms) 
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In this case, the superscripts denote the leading order 
terms (LO) while the subscript distinguishes the kind 
of interaction: spin-orbit (SO), spin(a)-spin(b) (SS) or 
spin(a)-spin(a) (S 2 ). We specified the initial spin of the 
particles by using a dimensionless parameter Xa € [0, 1] 
and the two spherical angles 9 a and <f> a - The initial spin 
of the particles is given by 

s o (0) = XaTO^cos^aSin^i + sin^sin^y + cos# a z), 

where x, y and z are the unitary basis vectors in Carte- 
sian coordinates. Using the Hamiltonian 0, the equa- 
tions of motion are 



(5) 
(6) 
(7) 



The first term in ([2| is the Hamiltonian for n-particles 
interacting under Newtonian gravity 
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The explicit form of the PN terms in ^ and the equation 
of motion for three compact objects can be found in [31] 
(see also [3H [49]). The spinning part terms ^ for n 
compact objects are given in [3T] , 

We will refer as Newtonian, 1 PN, 2 PN and 2.5 PN to 
the equations of motion derived from -H^" 1 ' H^p + Hp^, 

tt(0) Ml) rr(2) , „(0) „(1) „(2) „(2.5) roc „ of> 

h n +"pn+- h pn and -"n +- h pn+- h pn+- h pn j respec- 
tively. We denoted as radiative Newtonian and radiative 

1 PN the equations of motion derived from +i/p 2 ' J ' ) 
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B. Chaos indicators 

There are a number of methods for diagnosing chaos 
in a dynamical system: the Poincare surface (for system 
with less than 3 degrees of freedom), the Kolmogorov- 
Sinai entropy, the characteristic Lyapunov exponent and 
the fractal basin boundary method, among others (see 
e.g. [54TI56) ). For our analysis, we employed the frac- 
tal basin boundary method and the Lyapunov exponent. 
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In this section, we describe the implementation of both 
methods. 

The basin boundary method and the Lyapunov expo- 
nents are able to characterize different types of chaos. 
The basin boundary characterizes the sensitive depen- 
dence on initial conditions for nearby orbits with differ- 
ent attractors on the phase-spaceQ However, it does not 
provide information about the orbit itself. Alternatively, 
the Lyapunov indicator characterizes the regularity or 
chaos of an orbit in a way that is independent from the 
attractor. 



ln(JV ( 5 ii ( J 4))/ln(l/(5„) for a set of S n > l lcs . The typical 
behavior of the resulting data is well approximated by a 
straight line. A linear regression method is used to fit 
a linear function to the resulting data (where 5 n = nl res 
and n € {1, 2, ... , 12}). Therefore, the value of the fitted 
function for 5 = gives an estimate value of the dimen- 
sion. We tested the method by computing the dimension 
of non- fractal one-dimensional sets (i.e., a regular basin 
boundary) and Sierpinski's gasket and carpet. 

2. Lyapunov exponent 



1. Fractal basin boundary 

A dynamical system usually exhibits a transient be- 
havior followed by an asymptotic regime. The time- 
asymptotic behavior defines attracting sets in the phase- 
space which are called attractors. The attractor can be 
periodic, quasi-periodic or chaotic. The closure of the 
set of initial conditions, which has a particular attrac- 
tor, is called basin of attraction. Basin boundaries for a 
typical dynamical system can be either smooth or fractal 
[S"6l IS"?] . A system with fractal basin boundaries is sensi- 
tive to initial uncertainty. Given e > 0, we call the initial 
condition point P* unsafe if there is another initial con- 
dition point P inside a neighborhood of size e centered at 
P* which converges to a different attractor. The fraction 
/ of unsafe points is related to e by 



n-d f 



(9) 



where D is the dimension of the phase-space and d/ is 
the dimension of the basin boundary (see [56, 58J for a 
detailed discussion). It is useful to define an uncertainty 
exponent a = D — dy, which satisfies the relation < 
a < 1. The value of a quantifies the presence of chaos 
in a basin boundary, a rj 1 for regular boundaries and 
a w for a chaotic region. 

There are many ways of measuring the dimension of a 
fractal |57|, 159] . We employed the box- dimension which 
is defined as 



\n(N s (A)) 
dbox := hm - . . . ■ 
<5^o ln(l/d) 



(10) 



where A is a nonempty bounded set in the d^- 
dimensional Euclidean space and Ns(A) is the smallest 
number of sets in a 5-cover of A. A (5-cover of A is 
a collection of sets of diameter S whose union contains 
A. For a physical or numerical experiment, the fractal 
basin boundary has finite resolution. A graphical rep- 
resentation of the basin boundary is given by an image 
formed by pixels of size l Ies . In order to estimate the di- 
mension of the basin boundary, we compute the quotient 



1 An attractor is a set towards which a dynamical system asymp- 
totically approaches in the course of time evolution (see e.g. |54|). 



Lyapunov exponents are another way for characteriz- 
ing chaos in a dynamical system (see e.g. |54j ) . For an 
autonomous n-dimensional system 



X = F{X), 



(11) 



a given reference solution X* and a nearby solution X = 
X* + SX. While the evolution of the difference SX is 
given by the linearized equation 



5X = J(X») • 5X, 



(12) 



where J (A"*) 



djFi(X^) is the Jacobian matrix of F 



evaluated at the reference solution X*. The solution of 



( 12 1 is given by [60] 



6X{t) = c Jt • SX(0) 



(13) 



The Lyapunov exponents are defined by the eigenvalues 
Aj(t) of the distortion matrix A := c (J+J )* 



A, := lim — lnAi(i). 

t— >oo Zt 



(14) 



Sensitive dependence on initial conditions exists when 
Xi > 0. Therefore, in order to distinguish between reg- 
ular and chaotic orbits, it is sufficient to compute the 
principal Lyapunov exponent X p := max(A^). Calculat- 
ing X p directly is computationally expensive. However, 
by using the evolution of the difference, SX, it is possible 
to obtain an estimation of X v 



X* 



lim lim \Z(t,SX(0)), 



t-too 5X(0)-K) 



where, 



^(0)):4ln(g|), 



(15) 



(16) 



and 5X(t) = \SX\. We will refer as a Lyapunov in- 
dicator to A* in order to distinguish it from the prin- 
cipal Lyapunov exponent and the Lyapunov function, 
X*(t,SX(0)). The norm, |<5X|, plays an important role 
in the computation of the Lyapunov indicator. Partic- 
ularly, in the case of general relativity, additional gauge 
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effects can affect the estimation of the Lyapunov indica- 
tor [23j [53 [61] . Results from literature were applied to 
general relativistic dynamic of test particles and to bi- 
nary systems in the center of mass frame. Additionally, 
the estimation of A p is often computed using the two- 
particle method instead of the evolution of the difference 
(see e.g. [53]). However, we did not noticed significant 
differences in the computation of A* given the gauge ef- 
fects. We present the results in terms of the Cartesian 
length of SX. 

In practice, the Lyapunov indicator is computed ap- 
proximating the limits, using a sufficient small SX(0) 
and sufficient large integration time. In our simulations, 
we used |<DT(0)| = 10~ 9 . T he integration time depends 



on the system (see Sec. III). The evolution is performed 
by solving numerically the system of equations ([5])-([7|, 
therefore X — (x a ,p a , s a ) T is a vector with 27 compo- 
nents. The right hand side of (11) is the vector 
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and the equation of motion for the difference SX 
(Sx a ,Sp a ,Ss a ) is 
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where we define the auxiliary operator 
33/ 



6=i j=i \ ax b 




(21) 



III. SIMULATIONS AND RESULTS 
A. Numerical methods 

We solved the equations of motion numerically using 
the GNU Scientific Library (GSL) [53] and the Ollip- 
tic code infrastructure [32] • We generated the right 
hand side (RHS) of the equations of motion ([5|-([7]) with 
Mathematica 7.0 [64] . 

The simulations were done using the embedded Runge- 
Kutta Prince-Dormand (8,9) method provided by the 
GSL. We used a scaled adaptive step-size control (see 
[63j for details). In our simulations, the error control for 
the dynamical variables (x a ,p a , s a ) was set to 10 -11 and 
to 10~ 6 for the variation variables (Sx a ,Sp a ,Ss a ). We 
used as an indicator of accuracy the estimate of the local 
error provided by the GSL routines and the conserva- 
tion of the Hamiltonian for the simulations which do not 
included the radiation term. 



For our stability analysis, we solved the system of equa- 
tions ([5])-([7]) and (17)-(19) simultaneously. In order to 
perform the simulations in an efficient way, wc adapted 
the number of variables depending on the problem. The 
parameters to consider are the number of bodies nb, the 
magnitude of the spin s a and whether wc are solving a 
planar or non-planar problem. As we mentioned before, 
the RHS of Eqs. (|])-((7| were generated by a Mathe- 
matica script. For the Newtonian case and 1 PN cases, 
it is convenient to compute the analytical expression for 



the RHS of Eqs. (18)-(20). The analytical expression for 
the 2 PN and the spinning cases produces large compu- 
tational source files. The source file for the 2 PN is of 
25 megabytes while the inclusion of the spin generates 
files over a hundred megabytes. The compilation and 
optimization of large files is not practical, since compil- 
ers run out of memory even in the 16 GB of RAM server 
that we tried. Therefore, for the 2 PN order and the spin- 
ning case we computed the RHS of Eqs. ( 18 )-( 20 1 numer- 
ically. In order to estimate the errors in the numerical 
computation of the Jacobian, we compared the results 
obtained for the Newtonian and 1 PN cases employing 
both methods (see Sec. [Ill A 1[ ). The system ([l8|-([20]) is 
the Jacobian matrix of (|17p. The problem reduces to the 
computation of 



dFj(X t ) 
dXJ ' 



(22) 



A simple 2nd, 4th and even 6th order finite difference 
method with a fix step size h does not produce accurate 
solutions. The main issue is the nature of the compo- 
nents, for example, the position variables require a dif- 
ferent optimal step-size than the momentum or the spin 
variables. We implemented an adaptive method based 
on the GSL differentiation routine. The method pro- 
duces a solution with an estimated error below 10 -8 (see 



Sec. IIIAl) 



An important issue in the numerical integration of a 
three-body system arises when two of the bodies are very 
close to each other. In the case of adaptive step size 
methods, it is necessary to reduce the step size in or- 
der to resolve properly the orbits in the close interaction 
phase. A usual approach to deal with this problem is to 
perform a regularization of the equations of motion, see 
e.g., |65H69) and references therein. However, in our sim- 
ulations, we included a different criteria. Regardless of 
the equation of motion employed, wc monitored the ab- 
solute value of each conservative part of the Hamiltonian 
([I]) relative to the sum of the absolute values 
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#PnI + ^l(-ffSpin) , 



The simulations stop when the contribution of the first 
post-Newtonian correction is larger than 10%. We will 
refer to H^° > 10% as strong interaction. Empirical re- 
sults showed that, for a similar configuration, the resulted 
waveform exhibits a growth characteristic of the merger 
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phase [35] • The dynamics after a strong interaction may 
lead to the merger of two of the bodies or to the escape 
of the lighter body. In order to determine safely whether 
there is a merger or not, it is necessary to employ full 
numerical relativistic simulation. 

Additionally, we considered the Newtonian escape en- 
ergy of each body in respect of the other two components 



pa 

& csc - — 



Pa 

m„ 



m* 



(24) 



where fi a := m a m* a , m* := X^ a m 6 , %*a ■= 

( TO a) _1 Eb^a m bXb and p* := J2b^ a P>>- The simulation 
stops when E% sc > for one of the bodies and the size of 
the system S := E a l^<*l ^ s 100 times the initial size. Dur- 
ing a close encounter of two of the bodies it is common to 
have a positive value of E® sc for a short time. However, 
imposing the second condition we exclude those cases. 
We performed several tests to estimate the numerical er- 
rors. Here we summarize the results of these tests. 




3.5 

logioW 



(a) 



(b) 



J I I I l_ 



(c) 



1. Numerical tests 

Our first test was done with Lagrange's triangle solu- 
tion using Newtonian equations of motion. In Lagrange's 
solution each body is sitting in one corner of an equilat- 
eral triangle (see e.g. [TO]). We set the sides of such trian- 
gle to 7"i2 = ^23 = r 3i = 100, the mass ratio to 1:10:100, 
and the eccentricity to zero. Then each body follows a 
circular orbit (with different radii) around the center of 
mass. This configuration is unstable and after 14 periods 
one of the bodies escapes from the system. The unsta- 
ble property is reflected in the Lyapunov indicator A*. 
Fi ' 
for 



ig. flj(a) shows log 10 (A*) computed using two methods 
ir the prescription of the RHS of Eqs. (fl7l-Jl9l). The 



first method is the analytic expression (labeled J a ) and 
the second, the numerical solution (labeled J„). The two 
methods produce solutions that cannot be distinguish- 
able within the plot. The maximum relative difference is 
4 x 10~ 6 . The Lyapunov indicator exhibits the charac- 
teristic behavior of a chaotic system. After t = 10 4 , A*(£) 
oscillates around a constant value. Fig. [I] (b) shows the 
relative variation of the Hamiltonian 



AH := 



H(0) - H(t) 
H(0) ' 



(25) 



for the same solution. The conservation of the Hamilto- 
nian has a similar behavior either using J a or J„. The 
Hamiltonian exhibits a variation of around 10~ 12 before 
the evolution stops. 

Fig. [I] (c) illustrates the norm of the estimated 
errors for Lagrange's triangle solution. By looking at 
the errors, the difference between the solution produced 
by J„ and J a becomes clear. The numerical solution 
generates errors which are two orders of magnitude larger 
than the analytical Jacobian. However, even when using 
J a , the accuracy is good, with errors below 1.5 x 10 -8 . 



FIG. 1. Newtonian Lagrange's triangle solution. The up- 
per panel (a) shows the evolution of the Lyapunov indicator 
Xp(t,6X(0)), the middle panel (b) shows the conservation of 
the Hamiltonian AH and the lower panel (c) shows the L2 
norm of the estimated errors. In every panel, the solid line 
denotes the solution obtained using the analytical Jacobian 
and the dashed line, the result provided by the numerical 
computation of the Jacobian (see Sec. Ill A I. 



A second test was done for a stable system. Using the 
1 PN equations of motion, we computed the Lyapunov 
indicator A*, the relative variation of the Hamiltonian 
AH and the L2 norm of the estimated errors for the 
Henon's criss-cross solution [3H [7TJ [75] . We used initial 
parameters given by 

fi(0) = 1.07590A 2 i, pi(0) = 3~ 3 / 2 • C^OGA" 1 ^ 
f 2 (0) = -0.07095A 2 x, pa(0) = -3~ 3 / 2 ■ l^dmX^y, 
x 3 (0) = -1.00496A 2 £, p 3 (0) = 3~ 3 / 2 • 1.03678A _1 y, 

where x, y and z are the unitary basis vectors in Carte- 
sian coordinates, and A is a scaling factor (for our sim- 
ulation A = 10). Notice that for this test we used the 
parameters given in [73] with the scaling factor A, and 
doing a change of variables from initial velocity to ini- 
tial momentum. Therefore, we are not including post- 
Newtonian corrections to the initial parameters. As in 
the previous case, A* does not show significant differences 
when computed using J a or J„ (see Fig. [2]). The value 
of A* decreases monotonically almost like a straight line; 
this is the characteristic behavior of a regular solution. 
The maximum relative difference in A* using J a and J n is 
8 x 10 -6 . The Hamiltonian is conserved with a maximum 
variation of 5 x 10~ n . The L2 norm of the estimated er- 
rors shows noisy errors using J n while J a exhibits a flat 
trend. 
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FIG. 2. 1 PN Henon Criss-cross solution. The upper panel (a) 
shows the evolution of the Lyapunov indicator A* (t, <5X(0)), 
the middle panel (b) shows the conservation of the Hamilto- 
nian AH and the lower panel (c) shows the L2 norm of the 
estimated errors. In every panel, the solid line denotes the so- 
lution obtained using the analytical Jacobian and the dashed 
line, the result provided by the numerical computation of the 
Jacobian. 



A third test was done with a known chaotic binary 
system. We used a chaotic eccentric configuration with 
mass ratio 3:2 and initial parameters: 

f 12 (0) = 50x 

£i (0) = -p 2 (0) = 0.0061644?) + 0.003616z, 

s\(0) = -0.43765200x - 0.11469240y + 0.294789605, 

s 2 (0) = -0.02443680i + 0.10324000y - 0.01104320z. 

A similar system, with a different unit convention, was 
studied in [7J] (Sec. IV-C1). Following [7J], we solved 
the equation of motion using up to 2 PN correction for 
the orbital part and leading order in the spin. As was 
emphasized in [74] . in this configuration, the two com- 
pact bodies are rather close to each other. In order to 
obtain an accurate physical description, it is necessary 
to include, for the orbital part, corrections higher than 
2 PN and, for the spinning part, higher than the leading 
order. The 1 PN contribution to the Hamiltonian (22) 



during a close encounter of the bodies can reach 50%. At 
this point it is important to recall one of the conclusions 
of [74] . When using up to 2 PN Hamiltonian with leading 
order spin contributions, there is no evidence of chaotic 
systems for configurations which are physically valid to 
the PN order considered. Therefore, for comparison, this 
test was performed without using the stopping criteria 
H?°>10%, 

The result of the third test is presented in Fig. [3] The 



FIG. 3. Spinning binaries with an eccentric orbit. Lyapunov 
function X*(t,SX(0)) (a), conservation of the Hamiltonian 
AH (b) and the L2 norm of the estimated errors, for the 
system with initial parameters given by ( |26[ ) (solid line) and 
the same configuration except for s*i(0) = (dashed line). 



system was evolved using the parameters given in (|26j), 
and for reference, a non-chaotic orbit with the same ini- 
tial condition but with Si(0) = was used. Fig. [3} (a) 
shows the Lyapunov indicator for the chaotic configura- 
tion (solid line) and the arithmetic mean computed for 
t > 10 6 . The estimated value of the Lyapunov indicator 
is (A*) = 10~ 5 (dotted line) with a standard deviation of 
1.3 x 10 -6 , corresponding to a relative variation of 13%. 
In contrast, the regular solution shows that A* decreases 
almost monotonically after t = 10 4 . In both simulations, 
the Hamiltonian is numerically conserved below 10 -8 (see 
Fig. pjj-(b)), and the estimated errors are of order 10~ 6 
for the chaotic solution and 10~ 10 for the regular solu- 
tion. The main source for the errors, in the case of the 
chaotic solution, is the numerical calculation of the Ja- 
cobian matrix. 

We performed, additional tests already presented in 
Sec. III-A of [39] ■ A direct comparison with the results in 
[3"9"] does not showed significant differences in the results. 
For conciseness, we do not display the results of the tests. 



B. Stability of hierarchical systems 

Here, we consider the strong perturbation of a binary 
compact object system due to a third smaller compact 
object. As a basic configuration, we studied a Jacobian 
system with mass ratio 10:20:1. The inner binary sys- 
tem had initial apo-apsis ( i.e. maximum separation of a 
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TT 



FIG. 4. Hierarchical system (see [39]). Initial configuration 
of the inner and external binaries. The initial momentum of 
the third body is given by considering the external binary as 
a Newtonian binary. Shown are the osculating orbital planes 
Ilin and Ilext for inner and external binary orbits. The two 
planes are inclined by an angle t. 



Keplerian orbit) r b (0) = 130 and eccentricity e b (0) = 0. 
Although there are methods in post-Newtonian dynam- 
ics to specify the initial parameters of a binary system 
with a given eccentricity (see e.g. [74H79] ). in this study 
the inner binary is strongly-perturbed by a third body. 
Therefore, it is necessary to account for additional effects. 
For simplicity, we set the initial parameters considering 
only the Newtonian dynamics of a non-perturbed binary, 
where the eccentricity refers to the Newtonian case. In 
this approach, we view the third compact body and the 
center of mass of the inner binary as a new binary (we 
will refer to it as the external binary). The bodies start 
from a configuration where the initial radial vector r 12 is 
perpendicular to initial vector position x 3 of the exter- 
nal body (see Fig. [I]). We denote the inclination angle 
between the osculating orbital planes Ilj n and II GXt by 1 
(see Fig. |4]). To be more specific, the initial position and 
momentum of each body is given by 

Xi = ^r b x - m 3 r 3 (y cos i + zsini), 
■^2 = 7^r r bX - m 3 r 3 (y cos i + zsini), 



£3 = TOfeT"3 (y cos i + z sin i) , 

mi 



Pi 
P2 

P3 = P3X 



-PbV ~ ^T b P3X, 
PbV - ^P3X, 



(27) 



where m& := mi 
p 3 = m b m 3 ^ < /{^ 



m 2 , Pb 



mira 2 \/(l - 
Notice that 



e-b)/(m b r b ) and 
we set the mo 



e-3)/r 3 . 

mentum in terms of the eccentricity and the apo-apsis. 
Therefore, a higher eccentricity given to the external bi- 
nary implies a stronger perturbation of the inner binary. 

The main goal of this study was to characterize the 
stability and chaos of a hierarchical system as function of 
the initial apo-apsis ?'3, the eccentricity e 3 for the external 



binary and the inclination angle l between the osculating 
orbital planes. We explored the influence of the post- 
Newtonian corrections. 



1. Final state survey 

An analysis of the asymptotic behavior of the system 
as function of three parameters r 3 , e 3 and 1 was per- 
formed. For a given osculating angle t, we produced 
a map which is a subset of the space of configurations 
( r 3; e 3)- We define three possible outcomes which char- 
acterize the asymptotic behavior of the system: the es- 
cape of one of the bodies, a strong interaction E^' > 10% 
and a system which remains stable up to tf = 2 x 10 7 . 
We assign a color to each outcome in the map of ini- 
tial configurations (r3,e3): dark gray (color online) for 
an escape, white for a strong interaction and black for a 
stable configuration. Strictly speaking, the result is not 
a basin boundary map since the outcomes are not attrac- 
tors in the phase space and the parameters (?"3,e3) are 



not coordinates of the phase space. However, Eqs. (27) 



gives the link between (r 3 , e 3 ) and the coordinates of the 
phase space. On the other hand, a system where one of 
the bodies escapes will be attracted to a hyper-plane in 
the phase space with one of the spatial coordinates at in- 
finity; a strong interaction may become an escape or the 
merger of two of the bodies. A merger can be represented 
by a hyper-plane where two of the bodies share the same 
values of position and momentum. For a stable system, it 
is hard to define a single attractor because it can exhibit 
a different asymptotic behavior for t > tf. For simplicity, 
we refer to the boundary between the various color in the 
map as basin boundary since it can capture some of the 
features of a basin boundary. 

The choice of tf — 2 x 10 7 is motivated by the dy- 
namics of a similar Jacobian system presented in |39j . 
Including up to 2.5 PN terms, the inner binary arrives 
to the merger phase after i mgr w 2.7 x 10 7 . The assump- 
tion is that if a conservative system remains stable until 
that time, then the corresponding radiative system will 
remain stable until the merger of the inner binary. We 
selected a shorter time (i.e., tf < t mgi ) due to compu- 
tational costs. However, a set of numerical experiments 
for different values of tf suggest that a non-stable orbit 
is manifest before tf. 

We explored the parameter space r 3 , e 3 using 4 regions 
defined by 



Ri 
R 2 

R 3 



[200,500] x [0,1], 
[200,350] x [0,0.5], 
[275,350] x [0.25,0.5], 
[295,305] x [0.35,0.45], 



Ar 3 = 1, 
Ar 3 = 0.5, 
Ar 3 = 0.25, 
Ar 3 = 0.05, 



Ae 3 = 0.0025; 
Ae 3 = 0.00125; 
Ae 3 = 0.000625; 
Ae 3 = 0.00025. 



The domain of each region is the Cartesian product and 
its resolution is denoted by Ar3 , Ae3 . Each map was 
done computing 300 x 400 orbits except for region R 3 
that contains 200 x 400 orbits. Fig. [5] shows the result 
for a Newtonian simulation and osculating angle i = 0. 
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FIG. 5. Fractal basin boundary (Newtonian). Result of the 
asymptotic behavior of initial conditions for the parameters 
T3,e3 and i — for region Ro (a), Ri (b), R2 (c) and R3 
(d). The black points denote stable orbits, dark gray points 
(blue color online) are escapes of the lighter body and white 
points denote a strong interaction of two of the bodies. 



TABLE I. Newtonian orbits. Here db ox is the box dimen- 
sion of the basin boundaries presented in Fig. [5] the columns 
'Regular', 'Escape' and '-Hipn > 10%' denote the percentage 
of simulation that have the corresponding final state. 



Region 


dbox 


Regular 


Escape 


-fflPN > 10% 


Ro 


1.76 ±0.015 


8.2% 


32.4% 


59.4% 


Ri 


1.84 ±0.013 


1.0% 


41.3% 


57.7% 


R> 


1.77 ±0.015 


0.0% 


24.7% 


75.3% 


Rs 


1.74 ±0.007 


0.0% 


27.2% 


72.8% 



Notice that R3 c R 2 C Ri C Ro, i.e. each Ri is a magni- 
fication of the region Ri-i (in the figure, the sub-regions 
are indicated by rectangles). Fig. [5] shows the character- 
istic behavior of a fractal set; every magnification reveals 
a more complex structure from which it is possible to 
distinguish some degree of self-similarity. Table [I] sum- 
marizes the quantitative results for the basin boundary. 
Notice that the box-dimension reflects the fact that the 
sets are not perfectly self-similar. We choose R\ and R2 
to perform additional simulations. 

A comparison using Newtonian and 1 PN orbits for re- 
gions i?i and i?2 is presented in Fig. [6j There is a clear 
difference in the set of stable orbits. In the Newtonian 
case, there are two disjoint sets where the 1 PN system 
exhibits a single set (compare the black point in panels 
(a) and (b) of Fig. (6]) . The magnification R 2 shows dif- 
ferences in the escape and strong interaction zones. The 



FIG. 6. Fractal basin boundary (Newtonian and 1 PN). Com- 
parison of regions Ri and R2 for Newtonian and 1 PN orbits. 
Panels (a) and (b) show R\ for Newtonian and 1 PN respec- 
tively. Similarly, panels (c) and (d) show R2 for Newtonian 
and 1 PN, respectively. 

1 PN case exhibits slightly different substructures (see 
panels (c) and (d) of Fig. [6]). Our comparison includes a 
set of simulations for leading order spinning 1 PN parti- 
cles. We employed maximally spinning particles Xa = 1 
and 9 a = 0,(f) a = except by 6 1 = tt/2, 6> 2 = 2n/2 
(see Eq. Q). Therefore, the initial spin of the particles 
[configuration- a) is given by 

si = m\x, 

s 2 = —m 2 x, (28) 
s 3 = mp, 

Opposite to the non-spinning case, the orbits for our con- 
figuration of spinning particles are not coplanar. How- 
ever, the change in the final states is negligible, the re- 
sulting map is almost indistinguishable to the 1 PN case 
(for brevity, we do not display the basin boundaries which 
are similar to panels (b) and (d) of Fig. |6j. The quan- 
titative analysis of the six basin boundaries under com- 
parison is displayed in Table [TTJ The difference on the 
substructures between the Newtonian case and the 1 PN 
case is reflected by the box-dimension (smaller value for 
the 1 PN case). In the region R±, there is a small change 
in the distribution of 'Regular', 'Escape' and strong in- 
teractions between Newtonian and 1 PN. Nevertheless, in 
region R 2 there are noticeable differences in the percent- 
ages of 'Escape' and strong interactions. The similarity 
between the spinning and the non-spinning 1 PN case ap- 
pears for the dimension and the distribution of the three 
outcomes (with minor differences in the percentages). 
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TABLE II. Comparison of Newtonian, 1 PN, radiative 1 PN 
and leading order spinning 1 PN orbits. Here column PN 
takes value for Newtonian, 0+1 for 1 PN and 0+1+2.5 for 
radiative 1 PN. The column Sp takes value for non-spinning 
particles and 1 for the spinning ones. The meaning of the 
remaining columns is similar to TableJI] See the text for details 
about the parameters of the spinning case (Fig.|6]and[7]shows 
the non-spinning cases). 





PN 


Sp 


dbox 


Regular Escape -ffipN > 10% 










1.84 + 0.013 


1.0% 


41.3% 


57.7% 


Ri 


0+1 





1.81 + 0.014 


1.4% 


40.3% 


58.3% 




0+1 


1 


1.81 + 0.014 


1.2% 


40.7% 


58.1% 




0+1+2.5 





1.81 + 0.014 


1.7% 


40.0% 


58.3% 










1.77 + 0.015 


0.0% 


24.7% 


75.3% 




0+1 





1.74 + 0.015 


0.0% 


18.0% 


82.0% 




0+1 


1 


1.74 + 0.015 


0.0% 


18.1% 


81.9% 




200 220 240 2W 260 300 320 340 



FIG. 7. Basin boundary of region Ri. The dynamic includes 
(0+1+2.5) PN terms. 



In order to investigate the influence of the gravitational 
radiation, we performed an evolution which contains the 
2.5 PN terms. The inclusion of PN terms is computa- 
tionally expensive. Therefore, we used a radiative 1 PN 
Hamiltonian (i.e., we include the Newtonian, 1 PN and 
2.5 PN terms in ([2])). For similar configurations, the dif- 
ference in the dynamics between the full 2.5 PN system 
and the radiative 1 PN is relatively small [39]. Fig. [7] 
shows the fractal basin boundary and Table |TT] the corre- 
sponding quantitative measures. The radiative term does 
not seem to produce a significant change in respect of the 
inclusion of 1 PN terms (compare Fig. [7] with Fig.[6}(b)). 
Particularly, the box-dimension of the boundary is the 
same, and there is only a small difference in the distribu- 
tion of the percentages of escape and strong interactions. 

The osculating angle c has a strong influence in the re- 
gion of the initial parameters R\. For the Newtonian 
case, we produced a set of eight additional maps for 
region R\ where the osculating angle takes the values 
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FIG. 8. Fractal basin boundary of region R± as function of 
the osculating angle t (Newtonian) . The meaning of the color 
is the same as in previous plots (see e.g. Fig. [6jl . From the 
top to the bottom and from the left to the right i takes the 
values 0, tt/16, tt/8, 3tt/16, tt/4, 5tt/16, 3tt/8, 7tt/16, tt/2. 




FIG. 9. Fractal basin boundary of region R\ as function of 
the osculating angle l (1 PN). The meaning of the color is the 
same as in previous plots (see e.g. Fig. [6|. From the top to 
the bottom and from the left to the right l takes the values 
0, tt/16, tt/8, 3tt/16, tt/4, 5tt/16, 3tt/8, 7tt/16, t/ 2 - 
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7T/16 tt/8 37T/16 tt/4 5k/16 3n/8 ln/16 n/2 
l 

FIG. 10. Uncertainty exponent a as function of the osculat- 
ing angle i for the region R\. The upper panel (a) shows the 
result for the Newtonian case, the solid line is the fitting func- 
tion afl t (i) = ttNt+^N+CN sin^Nt- The lower panel (b) shows 
the corresponding result for the 1 PN simulations, where the 
solid line is a™ (t) = apNi + &pn + cpn sin^pNt. The fitting 



parameters are given in Eq. ( 29 1 



i = 7i7r/16, n € {1, 2, ... , 8}. The results, including the 
case i = 0, are presented in Fig. [8] An analogous set of 
simulations was produced using the 1 PN equation of mo- 
tion. Fig. [9] shows the corresponding basin boundaries. 
Table |III| shows the quantitative results for both sets. 
Notice that for i > ir/4 the percentage of stable points 
considerably increases for the 1 PN case in respect of 
the Newtonian. The box-dimension in both cases has a 
growing-oscillatory behavior. Fig. [10] shows the results 
for the uncertainty exponent a. A relatively simple func- 
tion cvfit(i) = ait + bi + CiSmipii fits the data. The fit 
parameters are 



a N = -0.059 ± 0.0088, a 1PN = -0.097 ± 0.0065, 

b N = 0.144 ± 0.0079, 6i PN = 0.183 ± 0.0059, 

c N = 0.026 ± 0.0062, c 1PN = 0.018 ± 0.0047, 

<p N = 7.1 ± 0.27, <^ipn =6.3 ±0.27. 



(29) 



The slope of the linear part of the 1 PN is 1.6 times 
larger than the Newtonian one. However, the oscillatory 
part is nearly 0.7 times the 1 PN value for the Newtonian 
case. This result shows that the chaotic properties of this 
hierarchical configuration increases in both cases, reach- 
ing the maximum at i — -k/2. From the basin bound- 
ary figures, we can notice an increasing number of unsafe 
points. There is an evident difference for i = n/2 between 
the Newtonian and 1 PN basin boundary (comparing the 
bottom-right panel of Figs. [8] and [9]). We analyze this 
particular case in the next section. 
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FIG. 11. Semi-major axis and Lyapunov function for initial 
parameters i — n/2,rz = 340, es = 0. The upper panel (a) 
shows the semi-major axis a of the inner and external binaries 
as function of time for the Newtonian case. The middle panel 

(b) is similar to (a) but for the 1 PN case. The lower panel 

(c) shows the evolution of the Lyapunov indicator for both 
cases. Notice that the quantities are given in logio scales. 



2. Analysis of orbits 

We used the basin boundaries as a guide for a detailed 
study of specific orbits. Given the region an osculat- 
ing angle t, two different basin boundaries n(r3,es) and 
^*( r 3j e 3): it is possible to produce a new map D that 
shows the differences between f2(r3,e3) and ft*(r^,e3). 
We constructed the map as follows; D(r 3 , e 3 ) takes value 
1 if 0(7-3,63) = il*(r^,el) in a neighborhood of (r 3 ,e 3 ) 
of size at least 6dr3 x 6de3, i.e., for — r% + idr 3 , 
63 = e3 + zde 3 and i G {~3, —2, . . . , 3}. If the previ- 
ous condition is not satisfied, then the D(r 3 ,es) is set to 
0. D{r^,e^,) provides a way for identifying zones where 
the two basin boundaries differed in a relatively large 
region. From the analysis of this map and the basin 
boundaries, we explored several initial parameters where 
we compared the orbits for different configurations of the 
PN equations of motion. In the following, we describe 5 
representative comparisons. 

For the region R\ and 1 = 7r/2 there is a clear difference 
between the Newtonian and 1 PN evolution (i.e., the dif- 
ference between the bottom-right panel of Figs. [8]and|9|. 
The bottom-right zone of parameters leads to a different 
outcome. In the Newtonian case, the simulation with pa- 
rameters e-i < (r 3 — 310)/320 ends when f/iPN > 10%. 
However, for the same set of parameters the 1 PN case 
remains stable. For the Newtonian case, there is a fast 
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TABLE III. Measures of the basin boundaries of Figs. [8] and [9] The osculating angle is denoted by l. The meaning of the 
remaining columns is similar to Table |ll| 





Newtonian 


1 PN 


6 


dbox 


Regular 


Escape 


H 1PN > 10% 


dbox 


Regular 


Escape 


H 1PN > 10% 





1.84 ±0.013 


1.0% 


41.3% 


57.7% 


1.81 ±0.014 


1.4% 


40.3% 


58.3% 


tt/16 


1.85 ±0.013 


1.0% 


43.0% 


56.0% 


1.82 ±0.014 


0.7% 


41.6% 


57.7% 


7r/8 


1.88 ±0.013 


1.1% 


49.2% 


49.7% 


1.84 ±0.013 


0.9% 


48.2% 


50.9% 


37r/16 


1.92 ±0.012 


1.9% 


62.3% 


35.8% 


1.90 ±0.012 


1.5% 


58.0% 


40.5% 


tt/4 


1.91 ±0.012 


2.7% 


67.9% 


29.4% 


1.91 ±0.012 


2.2% 


64.2% 


33.6% 


57r/16 


1.90 ±0.012 


2.7% 


74.3% 


23.0% 


1.91 ±0.012 


3.2% 


67.9% 


28.9% 


3tt/8 


1.89 ±0.013 


2.2% 


77.5% 


20.3% 


1.92 ±0.012 


4.3% 


70.8% 


24.9% 


7tt/16 


1.94 ±0.012 


1.5% 


73.9% 


24.6% 


1.94 ±0.012 


6.1% 


69.3% 


24.6% 


7r/2 


1.98 ±0.011 


0.1% 


57.1% 


42.8% 


1.97 ±0.011 


9.3% 


67.4% 


23.3% 



growth on the eccentricity of the inner binary induced by 
the external body (see Fig. [Tl]-(a)). This effect is pro- 
duced by a resonance between the orbital dynamic of the 
inner and external binary (see [SO] for a detailed descrip- 
tion). For the 1 PN case the orbits remain stable during 
the simulation (see Fig.[lTj(b)). In both cases, the evolu- 
tion of the Lyapunov indicator does not showed evidence 
of exponential growth behavior (see Fig. 11 -(c)). Notice 



a=1,b=2 



a=1,b=3 a=2,b=3 



that in the Newtonian case, the strong interaction is be- 
tween the components of the inner binary. The behavior 
of the Lyapunov function suggests that the strong inter- 
action between the heavy components of the system does 
not follow chaotic trajectories. 

For the region R\ and i = 0, there is a significant 
difference between the Newtonian and the 1 PN basin 
boundaries for stable orbits (i.e., the black dots in panels 
(a) and (b) of Fig. [6]). The Newtonian configuration re- 
mains stable without noticeable instability. On the other 
hand, the coupling between the orbits in the 1 PN case 
results in the escape of the third body (see panels (a) 
and (b) of Fig. 12 1. The Lyapunov function shows for 



1 PN dynamics the characteristic behavior of escape or- 
bits. During a close encounter between the lighter body 
and one of the heavy components of the inner binary, 
the difference vector has an exponential growth which 
is followed by a regular behavior. The regular behavior 
is expected for an uncoupled binary-single body system 
(see the kink in the Lyapunov function of Fig. [l2|(c)). 
The exponential growth for the difference vector is an 
indication of sensitivity to initial conditions. In general, 
the resulting escape orbit is chaotic; a small change in 
the initial condition produces a significant change in the 
escape direction. 

A comparison between 1 PN and 2 PN dynamics 
showed that almost identical orbits are produced by an 
evolution where the lighter body has a quick encounter 
with the inner binary which is followed by an escape or 
a strong interaction. Fig. [13] displays a comparison be- 
tween the Newtonian, 1 PN and 2 PN evolutions for 11 
nearby configurations. The initial parameters are t = 



o 



tu 
a 



3 


- Newtonian 


11,111,111! 


111,1 


2.7 








2.4 






2.1 














1.8 


" I I I I I 




■ ■ ■ 1 ■" 


4.8 


■ 1 PN ' 


1 1 1 1 1 1 1 1 1 1 1 




4.2 








3.6 






f ! 


3 








2.4 








1.8 


= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 J 1 1 1 1 1 — 



(b) 



Newtonian 



1 PN 




FIG. 12. Relative separation and Lyapunov function for ini- 
tial parameters i = 0, rz = 320, ez = 0. The upper panel (a) 
shows the relative separation r a b between particles as func- 
tion of time for the Newtonian case. The middle panel (b) 
is similar to (a) but for the 1 PN case. The lower panel (c) 
shows the evolution of the Lyapunov indicator for both cases. 



and e 3 = 0.1 with r 3 € {269.95, 269.96, . . . , 270.05}. The 
relative change in the initial separation r 3 between the 
configurations is ~ 0.004%. The Newtonian evolution 



(Fig. 13 (a)), exhibits a chaotic behavior. Every simula- 
tion ends with the escape of the lighter body. There are 
significant differences in the trajectories. On the other 
hand, the 1 PN and 2 PN evolutions (panels (b) and 



(c) of Fig. 13 ) produce a quick escape which are similar 
for each initial condition. The maximum final difference 
between the orbits is 1.3% for the 1 PN case and 1.5% 
for the 2 PN simulations. Moreover, the final difference 
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FIG. 13. Evolution of the coordinate position | a?3 1 for 11 
nearby configurations. Initial parameters l = 0, ez = 0.1 and 
r 3 G {269.95, 269.96, . . . , 270.05}. The upper panel (a) shows 
Newtonian case, the middle panel (b) shows the 1 PN evolu- 
tion and the lower panel (c) the 2 PN case. 



FIG. 14. Relative separation and Lyapunov function for ini- 
tial parameters i = 0, rz = 280, ez = 0. The upper panel (a) 
shows the relative separation r a b between particles as function 
of time for the 1 PN case. The middle panel (b) is similar 
to (a) but for the 2 PN case. The lower panel (c) shows the 
evolution of the Lyapunov indicator for both cases. 



between the 1 PN and 2 PN reference evolution r 3 = 270 
is 6%. 

However, for successive encounters the evolution is sig- 
nificantly different. Fig. 14 shows the resulted orbits for 
initial parameters t = 0,r 3 = 240, e 3 = 0. For these 
initial parameters, the resulting outcome is the escape 
of the lighter body. Nevertheless, in the 1 PN evolu- 
tion there is an ejectiorj^] of the lighter body before the 
final escape. In contrast, the ejection is not present in 
the 2 PN orbit (compare panels (a) and (b) of Fig. 14). 
There is a significant difference in the evolution of the 
Lyapunov function. For the 1 PN case, there is a dou- 
ble kink which corresponds to the close encounters before 
and after the ejection. Nevertheless, the 2 PN evolution 
presents a single long kink (see Fig. 14 -(c)). 

Fig. [15] shows a comparison between the Newtonian, 
1 PN and 2 PN evolutions for 11 nearby configurations. 
The initial parameters are l = and e 3 = with 
r 3 e {279.95, 279.96, . . . , 280.05}. Similarly to the case 
presented in Fig. [l3j the relative change in the initial 
separation r 3 between the configurations is ~ 0.004%. 
The set of configurations includes the solution presented 
in Fig. 
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r 3 = 280. The evolutions exhibit sensitive de- 
pendence on initial conditions. In the three cases, Newto- 
nian, 1 PN and 2 PN, relatively nearby initial parameters 



2 In the three-body problem literature, an ejection refers to a tem- 
poral large separation of one of the bodies which finally returns 



produce a significant change in the orbits. Moreover, the 
inclusion of PN corrections produced a noticeable change 
in the trajectories. 

For the parameters that we have studied, the inclusion 
of spin correction has a small influence in the final out- 
come. However, like in the 2 PN case, for some of the 
initial parameters, the evolution is completely different. 
For the initial parameters l = 0, r 3 = 320 and e 3 = 0.39, 
the 1 PN orbit ends with the strong interaction between 



particles 1 and 3 (see Fig. 16 (a)). In contrast, the spin- 
ning leading order 1 PN for the spin configuration- a ( |28[ ) 
results with the escape of particle 3 (see Fig. [l6[(bj). 
For the escaping orbit, the evolution produces a Lya- 
punov function which after t = 10 5 decreases linearly. 
The spinning case shows initially the same trend. How- 
ever, the strong interaction produces a short kink at the 
end of the evolution (see Fig. 16 -(c)). 

We explored a different set of spin parameters. Here 
we present one of the results of configuration-b given by 



s 2 = m^x, 
S3 



(30) 



m\i). 



Fig. [17] shows the result for the initial parameters 1 = 
7r/4, r 3 = 242, e 3 = 0.2 for the 1 PN and spinning leading 
order 1 PN cases. In this evolution, the non-spinning case 
results in the escape of particle 3. For the spinning parti- 
cles, the result is the strong interaction between particles 
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FIG. 15. Evolution of the coordinate position | a?3 1 for 11 
nearby configurations. Initial parameters l = 0, e$ = and 
r 3 G {279.95, 279.96, . . . , 280.05}. The upper panel (a) shows 
Newtonian case, the middle panel (b) shows the 1 PN evolu- 
tion and the lower panel (c) the 2 PN case. The bold-black 
line is the reference solution presented in Fig. |14| 



a=1,b=2 



a=1,b=3 a=2,b=3 




(a) 



FIG. 16. Relative separation and Lyapunov function for ini- 
tial parameters i = 0, rz = 320, ez — 0.39. The upper panel 
(a) shows the relative separation r ab between particles as 
function of time for the 1 PN case. The middle panel (b) 
is similar to (a) but for the leading order spinning 1 PN 
configuration- a (initial spin given by (281). The lower panel 
(c) shows the evolution of the Lyapunov indicator for both 
cases. 



FIG. 17. Relative separation and Lyapunov function for ini- 
tial parameters i — n/4,r3 = 242, e$ = 0.2. The upper 
panel (a) shows the relative separation r a b between parti- 
cles as function of time for the 1 PN case. The middle panel 
(b) is similar to (a) but for the spinning leading order 1 PN 
configuration-b (initial spin given by ( 30 1 ) . The lower panel 



(c) shows the evolution of the Lyapunov indicator for both 
cases. 



TABLE IV. Evolution of the eccentricity of the inner binary 
eb- Here we refer to the orbits in terms of the figure given 
in the column 'Fig.', e(,(0) is the initial eccentricity of the 
inner binary and et{tf), the final eccentricity computed before 
the corresponding simulation finish, (e^) and a(et,) are the 
arithmetic mean and standard deviation of et, respectively. 
The column 'FS' refers to the final state of the evolution. 





(b) 


# 


Fi 


y 


e b (0) 


e b (t f ) 


(e b ) 


o~(e b ) 


FS 






1 




11 


(a) 


0.004 


0.493 


0.112 


0.1453 


tflPN > 10% 






2 




11 


(b) 


0.040 


0.040 


0.039 


0.0003 


Regular 




3 




12 


(a) 


0.021 


0.030 


0.025 


0.0047 


Regular 






4 




12 


(b) 


0.036 


0.281 


0.252 


0.0556 


Escape 






5 




14 


(a) 


0.042 


0.140 


0.117 


0.0225 


Escape 




(c) 


6 




14 


(b) 


0.062 


0.104 


0.104 


0.0074 


Escape 






7 




16 


(a) 


0.030 


0.155 


0.185 


0.0332 


tflPN > 10% 






8 




16 


(b) 


0.137 


0.184 


0.184 


0.0049 


Escape 




9 




17 


(a) 


0.095 


0.237 


0.225 


0.0343 


Escape 




10 




17 


(b) 


0.024 


0.004 


0.071 


0.0299 


H 1PN > 10% 



1 and 3. The Lyapunov function shows the characteris- 
tic exponential growth of the difference vector during the 
close interaction. 

Tabic |IV| shows the evolution of the eccentricity of the 
inner binary ej,. The eccentricity is computed by calculat- 
ing numerically the value of the apo-apsis and peri-apsis 
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of the inner binary i.e., computing a sequence of local 
maxima and minima of r\i- The eccentricity is given by 



Notice that the initial value of eb in every case is relatively 
small but not zero due to the perturbation of the third 
body. However, the final value of the eccentricity depends 
on the final state of the evolution. For regular orbits, 
the change is relatively small (see the standard deviation 
of rows 2 and 3 in Table |Tv| . On the other hand, the 
final eccentricity of the escaping orbits is relatively large 
during the whole evolution (compare the final eccentricity 
and the arithmetic mean of the corresponding rows). The 
orbit with resonance exhibits the largest increment in 
eccentricity. The final value before the evolution stops is 
close to 0.5. For PN evolutions, we did not find evolutions 
with this kind of resonance. However, a different set of 
parameters may lead to similar behavior. For a strong 
interaction between two of the bodies, the resulting final 
eccentricity is not meaningful since there are two possible 
outcomes in the full evolution scenario. One possibility 
is an ejection or escape of the lighter body leading to a 
increment in the eccentricity of the inner binary similar to 
the cases presented. The other possibility is a successive 
merger of the three bodies. 

IV. DISCUSSION 

We performed a numerical study of stability and chaos 
of a hierarchical three compact object configuration. The 
configuration is composed by a binary system which is 
perturbed by a third, lighter body positioned initially 
far away from the binary. For a given inner binary, we 
explored the influence of the third body taking as free pa- 
rameters the angle of the osculating orbital planes, the 
initial apo-apsis and eccentricity of the third body. Using 
the basin boundary method, a total of 2,960,000 simula- 
tions were analyzed. From the total, 47.3% corresponded 
to Newtonian evolutions, 40.54% to (0+1) PN simula- 
tions, 4.05% to (0+1+2.5) PN evolutions and 8.11% are 
simulations of spinning particles including leading order 
in spin and 1 PN terms in the orbital part. 

A total of twenty-five basin boundary maps were pro- 
duced. The basin boundaries exhibit the characteristics 
of fractal sets: some degree of self-similarity and an in- 
creasing complexity under magnification. Fractal basin 
boundaries are an indicator of chaos in dynamical sys- 
tems. By measuring the fractal dimension of the basin 
boundaries it is possible to determine the uncertainty 
exponent a. The property of the exponent is that for 
a « the dynamical system is chaotic and for a = 1 
the system is regular. The values of a for the fractal 
basin boundaries considered here are between 0.02 and 
0.26. The 1 PN set of data produces an exponent slightly 
larger than the Newtonian case (the maximum difference 
is 1.6%). On the other hand, the osculating angle i has 



a strong influence in the uncertainty exponent. The ex- 
ponent decreases as a linear oscillatory function in the 
range [0,7r/2]. The difference between the planar case 
(i = 0) and the case where the third body starts from a 
direction perpendicular to the orbital plane of the inner 
binary (i = tt/2) is 7.1% for the Newtonian simulations 
and 8.1% for the 1 PN case. 

In addition to the uncertainty exponent, we quantified 
the percentage of stable orbits, escapes and strong inter- 
actions. The configurations selected contain between 0% 
and 9.3% of stable orbits. The distribution of escapes 
are between 40.3% and 77.6%. The strong interactions 
oscillate between 20.3% and 82%. A remarkable case 
was found at l — tt/2 where the Newtonian simulation 
presents a resonance which drops the number of stable or- 
bits to 0.1%. Opposite to the Newtonian dynamic 
including 1 PN or higher corrections eliminates the res- 
onance. The number of stable orbits increases to 9.3%. 
The presence of resonances in the three-body problem is 
fundamental for understanding the chaotic properties of 
the system, see e.g. [T51 [STJ HO] and references therein. 

By looking at specific orbits, it is possible to notice 
some of the different couplings between the inner and 
external binaries. For most of the orbits, the 1 PN terms 
have the strongest effect on the dynamics. The inclusion 
of the spinning particles or 2 PN corrections produces 
small differences in the orbits when the lighter body has 
a quick encounter with the inner binary which is followed 
by an escape or a strong interaction. In that cases the 
orbits are essentially identical. However, the coupling of 
the orbits for some of the configurations can produce a 
significant change in the final outcome. We presented 
five representative examples. 

The spin can produce a significant change in the orbits 
producing precession of the orbital planes. However, the 
change on the dynamics is not strong enough to modify 
the final outcome. On the other hand, gravitational radi- 
ation, in general, produces a small change in the energy 
to be significant in the short time scale of the evolution. 
Similarly to the spinning case, gravitational radiation has 
a small effect in the asymptotic behavior. 

Considering the scenario of a binary system in a galac- 
tic core or a region with high density of compact objects, 
we expect to have the following results. Between a 40% 
and 70% probability that the lighter body escapes the 
system depending on the relative orbital planes. The re- 
manent binary may result in eccentricities between 0.1 
and 0.2. Between 23% and 58% of the bodies may lead 
to a merger with one of the components of the inner bi- 
nary; contributing to the growth of the compact objects. 
The probability of finding a body in a stable orbit is 
between 1% and 10%, producing a perturbation which 
leads to small eccentricities of the inner binary of around 
0.03. Nevertheless, the combined effect of many bodies 
of comparable size can significantly change the results. 

More general statements about the chaotic behavior 
of three compact bodies are possible but require an ex- 
tensive parameter study. Other configurations include, 
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for example, a characterization of the stability and chaos 
based on the mass ratios, inner binary separation, magni- 
tude and direction of the spin or initial eccentricity of the 
inner binary. We consider this a topic for future study. 
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